home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / mathpack / dgesl.f < prev    next >
Text File  |  1989-08-17  |  3KB  |  118 lines
  1.       subroutine dgesl(a,lda,n,ipvt,b,job)                              
  2.       integer lda,n,ipvt(1),job
  3.       double precision a(lda,1),b(1)
  4. c
  5. c     dgesl solves the double precision system
  6. c     a * x = b  or  trans(a) * x = b
  7. c     using the factors computed by dgeco or dgefa.
  8. c
  9. c     on entry
  10. c
  11. c        a       double precision(lda, n)
  12. c                the output from dgeco or dgefa.
  13. c
  14. c        lda     integer
  15. c                the leading dimension of the array  a .
  16. c
  17. c        n       integer
  18. c                the order of the matrix  a .
  19. c
  20. c        ipvt    integer(n)
  21. c                the pivot vector from dgeco or dgefa.
  22. c
  23. c        b       double precision(n)
  24. c                the right hand side vector.
  25. c
  26. c        job     integer
  27. c                = 0         to solve  a*x = b ,
  28. c                = nonzero   to solve  trans(a)*x = b  where
  29. c                            trans(a)  is the transpose.
  30. c
  31. c     on return
  32. c
  33. c        b       the solution vector  x .
  34. c
  35. c     error condition
  36. c
  37. c        a division by zero will occur if the input factor contains a
  38. c        zero on the diagonal.  technically this indicates singularity
  39. c        but it is often caused by improper arguments or improper
  40. c        setting of lda .  it will not occur if the subroutines are
  41. c        called correctly and if dgeco has set rcond .gt. 0.0
  42. c        or dgefa has set info .eq. 0 .
  43. c
  44. c     to compute  inverse(a) * c  where  c  is a matrix
  45. c     with  p  columns
  46. c           call dgeco(a,lda,n,ipvt,rcond,z)
  47. c           if (rcond is too small) go to ...
  48. c           do 10 j = 1, p
  49. c              call dgesl(a,lda,n,ipvt,c(1,j),0)
  50. c        10 continue
  51. c
  52. c     linpack. this version dated 08/14/78 .
  53. c     cleve moler, university of new mexico, argonne national lab.
  54. c
  55. c     subroutines and functions
  56. c
  57. c     blas daxpy,ddot
  58. c
  59. c     internal variables
  60. c
  61.       double precision ddot,t
  62.       integer k,kb,l,nm1
  63. c
  64.       nm1 = n - 1
  65.       if (job .ne. 0) go to 50
  66. c
  67. c        job = 0 , solve  a * x = b
  68. c        first solve  l*y = b
  69. c
  70.          if (nm1 .lt. 1) go to 30
  71.          do 20 k = 1, nm1
  72.             l = ipvt(k)
  73.             t = b(l)
  74.             if (l .eq. k) go to 10
  75.                b(l) = b(k)
  76.                b(k) = t
  77.    10       continue
  78.             call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
  79.    20    continue
  80.    30    continue
  81. c
  82. c        now solve  u*x = y
  83. c
  84.          do 40 kb = 1, n
  85.             k = n + 1 - kb
  86.             b(k) = b(k)/a(k,k)
  87.             t = -b(k)
  88.             call daxpy(k-1,t,a(1,k),1,b(1),1)
  89.    40    continue
  90.       go to 100
  91.    50 continue
  92. c
  93. c        job = nonzero, solve  trans(a) * x = b
  94. c        first solve  trans(u)*y = b
  95. c
  96.          do 60 k = 1, n
  97.             t = ddot(k-1,a(1,k),1,b(1),1)
  98.             b(k) = (b(k) - t)/a(k,k)
  99.    60    continue
  100. c
  101. c        now solve trans(l)*x = y
  102. c
  103.          if (nm1 .lt. 1) go to 90
  104.          do 80 kb = 1, nm1
  105.             k = n - kb
  106.             b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
  107.             l = ipvt(k)
  108.             if (l .eq. k) go to 70
  109.                t = b(l)
  110.                b(l) = b(k)
  111.                b(k) = t
  112.    70       continue
  113.    80    continue
  114.    90    continue
  115.   100 continue
  116.       return
  117.       end
  118.